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Abstract 

A vector field splitting approach is discussed for the systematic derivation of numerical propagators for 
deterministic dynamics. Based on the formalism, a class of numerical integrators for Langevin dynamics are 
presented for single and multiple timestep algorithms. 
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1 Introduction 



The design of efficient and stable propagators is a central issue in the numerical study of classical and quantum 
systems. The way the state is propagated in time affects the quality of the simulated trajectory and static and 
dynamical ensemble averages, as much as the capacity to efficiently sample the available phase space. 

Molecular Dynamics (MD) is the reference technique for the study of condensed matter systems [lj. Al- 
though, in principle, time discretization introduces an error which affects the statistics of the simulated systems, 
from the point of view of numerical control, MD offers a number of advantages. The underlying conservative 
and time-reversible nature of the Hamiltonian equations allows to monitor the quality of the dynamics as a 
function of the time step and the implemented interaction potential. Despite energy conservation is never at- 
tained in practice, keeping the energy fluctuations small is a necessary, although not sufficient, condition to 
avoid statistical bias in the computed averages. 

The lack of drifts in the energy has long been recognized as a key requirement to control long-time stability. 
In mechanical terms, long-time stability follows from the symplectic nature of the propagator, a distinguishing 
feature of Hamiltonian dynamics [2]. In statistical terms, symplecticity implies exact conservation of measure 
in phase space, i.e. the possibility of applying Gibbs statistical mechanics to a well-defined and conserved 
ensemble of physical states. Moreover, symplecticity and time-reversibility are the basis to combine MD with 
Monte Carlo, in the so-called Hybrid Monte Carlo method [3], where the acceptance test removes any systematic 
bias due to the usage of a large timestep. Another important benefit of conservative dynamics is the possibility 
to employ ad-hoc quasi-Hamiltonian dynamics to sample ensembles different from the microcanonical, a popular 
choice being the Nose-Hoover dynamics [TJ by retaining the same pleasant features of the energy conserving 
dynamics. 

One of the most popular numerical propagators originates from the pioneering work of Verlet (4| who 
employed an algorithm first introduced by Stormer [5] . The distinguishing features of the Verlet (also known as 
Stormer- Verlet) algorithm are its simplicity, as compared to predictor-corrector or Runge-Kutta methods [2], 
time-reversibility and the lack of numerical drifts. Therefore, in spite of its limited accuracy, being quadratic in 
the timestep, the Verlet propagator stands as the reference algorithm for MD. Some years after its discovery, a 
systematic derivation based on an operatorial splitting approach was proposed [6]. By considering the case of a 
separable Hamiltonian, the so-called Trotter factorization [7] was applied to derive the Verlet algorithm and its 
equivalents, in particular the Velocity Verlet (VV) version. A crucial benefit of such operatorial approach is its 
generalization to a multiple time step (MTS) scheme [5], which extends the single time step (STS) one. In the 
MTS, the advance in time of the phase space state follows the temporal evolution of forces in an asyncronous 
way, i.e. by treating fast and slow interactions at different levels. The efficiency of simulation can thus be 
improved up to one order of magnitude in favourable circumstances. 

The operatorial approach is very elegant and fruitful. However, it can be cumbersome in the systematic 
derivation of propagators in different contexts. A typical case is provided by non-separable Hamiltonians, finding 
wide application in the description in generalized coordinates or in presence of velocity-dependent forces (see 
e.g. |9]). In this case, symplectic and robust numerical schemes, such as the Generalized Leapfrog, have been 
known for some time (2]. Another wide class of dynamics is represented by stochastic equations of motion, 
such as the all-famous Langevin equation. In this context, the operatorial route has again been applied by 
considering the stochastic noise as a time-dependent perturbation [10] or by evolving the state according to a 
Fokker-Planck propagator jlll [T2] . It has to be mentioned that for a first-order stochastic differential equation, 
a fourth-order accurate scheme has been derived jTT] based on a Runge-Kutta propagation of the deterministic 
component. However, such scheme requires to compute high order derivatives, thus being rather complicated 
for practical applications in condensed matter, and is not specifically designed to reduce to symplectic in the 
limit of zero friction. In general, due to the presence of both stochastic and deterministic forces, conventional 
operatorial calculus should be applied with some care[13 . As a matter of fact, a previous operatorial approach 
was shown to overestimate the accuracy of the numerical trajectory |10l 114] . 

Stochastic dynamics applied to condensed matter systems has recently received renewed attention for sev- 
eral reasons: i) recent emphasis on multi-scale methods is based on the simultaneous evolution of atoms and 
surrounding hydrodynamic fields, the latter been treated via Lattice Boltzmann dynamics. In this technique 
noise serves to effectively enforce fluctuation-dissipation relations in otherwise decoupled systems, thus requiring 
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stable and accurate numerical integrators |1 51 H6] ; ii) improving the stability of deterministic MTS algorithms 
to large timesteps can be achieved by attaching to each degree of freedom massive thermostats to overdamp the 
atomic motion [T7] or by stabilizing the motion by an overdamped Langevin friction. The notion of Langevin 
stabilization is not new, and has already been proposed in the literature in different forms |18U19] : iii) Langevin 
dynamics is an effective mean of attaching the system to a thermal reservoir, thus allowing to sample the 
canonical ensemble without the introduction of artificial dynamics with memory (e.g. Nose-Hoover). 

In the present paper we will consider a deterministic class of propagators, the Verlet being a particular 
case for separable Hamiltonian dynamics, derived from a different perspective, based on splitting the phase 
space vector field rather than by approximating the Liouvillean. The trajectory representation employed here is 
clearly equivalent to the operatorial one. However, in the operatorial approach one evaluates the involved time 
integrals at current times (generating so-called "forward" update) while in the trajectory representation one can 
evaluate such integrals with more general interpolation schemes. The present approach allows to establish a clear 
connection between the algorithm and eventual approximations supplementing the propagation of the trajectory. 
In the same spirit, the technique is extended to Langevin dynamics for both STS and MTS approaches. We 
will derive numerical integrators for Langevin dynamics accurate to second-order in the timestep and reducing 
to symplectic Verlet ones in the limit of zero friction. Our derivation treats momenta as natural elements of the 
algorithm, an aspect which has recently been debated [20]. Some numerical tests demonstrate the statistical 
consistency of the schemes for STS and MTS propagations. Moreover, the stochastic MTS scheme applied to a 
realistic system made of water molecules, therefore comprising intramolecular, excluded volume and electrostatic 
forces, is shown to achieve improved numerical stability and correct configurational statistics. 



2 Deterministic dynamics 

Let us consider the 67V phase space point x = {qi,Pi}i=i.3N associated to the autonomous equations of motion 
written as 

x = SV(x) (1) 

where S is a generic matrix with constant elements and V(x) a generic vector field. In case of Hamiltonian 
dynamics the vector field reads 

SV(x) = SdH(x) (2) 
where d = {^r, g^:}- <5 is the so-called symplectic matrix 

where the 3N x 3N block matrices and Id are the null and identity matrix, respectively. Therefore, S has 
a trivial inverse S^ 1 = S T = —S, expressing the fact that the symplectic transformation is norm preserving 
(||<Sff|| = HsID- For separable Hamiltonians V(x) reduces to the familiar form 

**(.)-(-$>) M 

where F and p refer to the 3N components of forces and momenta and m are the masses, assumed to be equal 
for the sake of simplicity. 

In order to integrate numerically eq.([lj, we consider the generic decomposition 

S = U + L (5) 

where, for Hamiltonian flow, S reduces to S and U and L specialize to 

MooO '-(-IS) <°> 
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We employ now a splitting technique to integrate the equations of motion. The idea is to approximate 
the vector field (JTJ as the superposition of two linearly independent vector fields, each one associated to the 
equations of motion 

fx = LV ^ < 7 > 

and 

% = uv ^ < 8 > 

Symbolically, each evolution operator is briefly indicated as ^ and ^ [21], where ^ refers to the complete 
evolution ([TJ , associated to the updates 



e h *x (9) 



and 



x ^e h ^x (10) 

together with the composite evolutions x — > e h ^ e h ^x a and x a — > e /l 3xe /l3 f r x a . The two operators ^ and j- 
do not commute in general, i.e. the so-called Lie bracket j^} ^ 0, implying that the distance between the 
evoluted points has error of order h 2 , i.e. e h ^e h ^x — (e h *>*e h tt + h 2 [-^, jj;])x + 0(h 3 ) [21;. On the other 

hand, by considering the symmetric composition ei a* ei a* this approximates the true trajectory to 
second-order accuracy in a timestep h. 

Associated with the two evolutions ( 9|10 1 we now consider either the exact solution of the dynamical equa- 
tions for each vector field, <f>L,h(x ) and <fiu,h(x ), or the maps based on the four elementary approximations to 
the time integrals L V(t)dt and U Jq V(t)dt as 

®L,h{xo) = x a + hLV 

®Lji{x ) = x a + hLV h 

$u,h( x o) = %o + hUV Q 

§u, h {x ) = x + hUV h (11) 

where the bar denotes evaluation of the integral at the upper extremum of the interval (requiring in principle the 
inversion of an implicit equation). Conversely, the operatorial route is always based on evaluations at the lower 
extremum, thus allowing for more flexibility in the design of numerical algorithms in the present treatment. 
Nevertheless, the vector field splitting and the approximated integrals are two independent sources of numerical 
error. 

The Euler-A propagator is defined by the composition 

X h /2 = ®U,h/2 ° ^L,h/2{x D ) (12) 

and the Euler-B scheme is defined by the composition 

%h/2 = &LM/2 ° ®U,h/2{Xo) (13) 

It is easily shown that Euler-B is the adjoint of Euler-A, where the adjoint of the map Mh is defined as 
M£ = Mz\ 0. A symmetrized algorithm can be constructed as 

Xh = $L,h/2 ° ®U,h/2 ° ®U,h/2 ° ®L,h/2(Xo) (14) 

or by interchanging U <-> L. Given the presence of pairs of the type o <& L and o in the consecutive 
applications of eq.(14l, the global propagator is based on evaluation of the integrals via the midpoint method, 
being accurate to quadratic order in the timestep. The scheme is symmetric with respect to time-inversion, 
since 

Xo = ®L-h/2 ° ®U,-h/2 ° ®U-h/2 o ®L,-h/2{Xh) (15) 
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For non-separable Hamiltonians, the propagator (14 1 reduces to the Generalized Leapfrog [2]. The updating 
scheme reads 



p* 


= Po 


hdH, 






= q 


h 

+ 2 






Ph 


= p* 


hdH . 





where p* indicates intermediate values of the momenta. 

For a separable Hamiltonian each update of the map Xh — 
resulting scheme given by the celebrated Velocity Verlet (VV) 



(16) 



l) c,h/2 4>u,h ° 4>cm/2(x ) becomes explicit, the 



P 

% 

Ph 



Po 



q a + h 



: F{qh) 



(17) 



The so-called Position Verlet propagator (PV) is similarly obtained by composing the maps as Xh — <foi,h/2 ° 

4>CM ° 4>U,h/2{x ). 

A central property of the Generalized Leapfrog and Verlet-based propagators is symplecticity. Generally 
speaking, the symplectic property applies to Hamiltonian dynamics, which obeys the relation |2 



j t s~ x j = s- 1 



(18) 



where J is the Jacobian of the transformation x(0) — > x(t). This property is related to the fact that the 
dynamics is invariant under canonical transformations. Conservation of measure follows directly, since eq.(18l 
implies |j J|| 2 = 1. 

An explicit calculation proves that the Euler-A scheme for Hamiltonian flow is symplectic. In fact, 



J 



OX 



h/2 

2 dx Q 



= i + hsn h/2 (i - hcn,^)- 1 = (i + hcn h/2 ){i - mn h/2 y l 



(19) 



where = Cl T = d 2 H and having using the fact that x a — x h / 2 — hCdH h / 2: so that a ^ /2 = 
(I - hCn h/2 )-\ Finally, 

JTS^J-S- 1 = 



I - 

p- 



-Sfi h / 2 (i - -Cfi h / 2 y 



s- 



—sQh/2(i — -C£i h / 2 y 



-s- 



-LQ, h / 2 ) 



-liT 



h/2- t -"/i/2 



-Q, h / 2 UQ, h / 2 



h 



£lh/2S ^ 



h/2 



x(i--cn h/2 )- 1 = o 



(20) 



A central property of symplectic propagators regards the presence of a shadow Hamiltonian exactly conserved 
by the numerical flow and written as a power expansion in the timestep, H(x) — H(x) + ^2 n h n G n {x)^ such 
that the associated vector field is equal to Air = x^ — xq = hSdH. An explicit expression for the leading G n (x) 
terms can be derived with leading terms hG 1 (x) and —hG 1 (x) for symplectic Euler-A and Euler-B schemes 
respectively. As a consequence, VV scheme has leading term of h 2 order [2]. Clearly, on the H — const 
manifold, measure is preserved, since d ■ Ax — 0, and one can define a shadow distribution function f(y, t) such 
that a Liouville equation is exactly satisfied 



f(y,h)-f(y,0) + ~Axdf(y,0)=0 



(21) 
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3 Langevin dynamics 



Let us consider the Langevin dynamics 

(22) 



q \ _ / p/m 
P J V F(q)-jp+ V (t) 



where rj(t) is a white noise, with 



(»?(*)) - o 

fa(f)»j(0} = 5 2 <5(* - (23) 

7 is the friction coefficient, g 2 = 2-fmkT, and the deterministic part arises from a separable Hamiltonian. The 
present treatment, applied to the case of a scalar constant friction, can be rapidly extended to the case of a 
position-dependent or tensorial friction, as used in the presence of hydrodynamics or confined diffusion. The 
equations of motion are written in compact notation as 

x = [S + mG)dH(x) + fj(t) (24) 

where 

G -<tl) *>=U) 

Different algorithms can be employed to integrate the Langevin dynamics, the most popular being the Euler 
and the Heun ones [22] whereas, in presence of interatomic forces, different algorithms have been proposed in 



the past j23( [Ml l25l l26l [27] . The numerical solution of (22 I presents, in general, non trivial mixing between 
deterministic and stochastic terms. This is shown by a standard argument [28 . Let us consider an elementary 
step written in Euler form, Sxt = f t + Wt, where f(x) = (S + mG)dH(x) and Wt = J fj(s)ds is the Wiener 
increment [30]. The Euler step is exact at sufficiently small time t and can now be nested in the integral of 



(24 1, leading to the formal solution 

f h 1 
x h ~x = J dt(f + df 5x t + -d 2 f Sx 2 + .... + rj t ) = 

= foh + W h + df foY + dfoJ o dtW t + 0{h b ' 2 ) (26) 

This expression shows that, in order to obtain the desired accuracy, a second order integrator for the determin- 
istic component is needed together with a proper handling of the stochastic forces. In fact, to second order, 
deterministic and stochastic terms mix up via a linear functional of the noise, dfo J Q dtWt, scaling as ft, 3 / 2 . 
Due to the separable Hamiltonian, direct inspection reveals that this term reduces simply to 1 jm J Q dtWt for 
positions and —7 dtWt for momenta, i.e. arising from inertial and Langevin forces only. In principle, the 

purely deterministic terms (foh + dfofo \ + ■■■), accounting for the effect of the operator e h ^ 9 (xo) 1 cannot be 
separated from the stochastic ones via the simplistic operatorial splitting, at the expense of missing the mixed 



term in eq.(26). However, it should be possible to decompose the vector fields by separating the evolution of 



the configurational degrees of freedom from that of the momenta, by considering the splitting 

x = (C + mG)dH{x) + r)(t) (27) 

and 

x = UdH(x) (28) 



Eq. (27 1 has solution given by the exact map 



MXo) S {e-^p 0+ ^F + Qh ) (29) 
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where the following stochastic integral appears 



Qh = e^ h / dse^s) = N(mkT{\ - e^' 1 )) 



(30) 



with Af(a 2 ) being a gaussian variable with zero mean and variance a 2 . 
We define the stochastic Euler-A as 

and the stochastic Euler-B as 

%h/2 = 4>CJi/2 ° 4>U : h/2(Xo) 

By symmetric composition, the Stochastic Velocity Verlet (SVV) is constructed as 

Xh = 1pC,h/2°4 l U,h°4 , C,h/2(x ) 

where we have collapsed 4>u,h/2 ° 4 l u,h/2 — <hi.h- The complete updating scheme reads 
P* = e-~< h / 2 p + (1 ~ 6 ^'^ F(q ) + ffW (mkT(l - e~^)) 



(31) 



(32) 



(33) 



qh = q a -\ — p 

m 



Ph 



(34) 



where Af^ and Af^ are two independent realizations of the process (301 and thus, the algorithm requires two 
extractions of random numbers per timestep. 



By reverse composition, the stochastic position Verlet (SPV) is constructed, reading 

Xh = 4>U,h/2 ° ipc,h ° <i>u,h/2{x ) 



(35) 



where it is easily shown that ipc.h/2 ° 4>c,h/2 — ip£,h, having collapsed the sum of independent gaussian terms 
into a single gaussian extraction. In explicit terms, SPV reads 



q 

Ph 

qh 



qo + 

2m 

e-~t h p + (1 ~ 6 F(q*)+Ar(mkT(l - e" 2 ^)) 



q + ^—Ph 

2m 



(36) 



The accuracy of the SVV and SPV algorithms is evaluated below. At first, let us consider the explicit 



solution of eq.(22 1 



x, 



q + J* Me'*' A [p + dt"e*"F( q (t"))} + i J* df f*' dH' e*" n {H>) 
e'^po + e-T* /* dt'e*'F(q{t 1 )) + J* dt 1 e*" r,{t') 

+ (^) Po+^ Jo dt'F(q(t')) (l - e-tf-O) + Q t ) \ 



e-*p + /„* dt'e*'F(q{t')) + Q t 



(37) 



where we used the rule J Q dt' J Q dt" = dt" f*„ dt' . The variable Q t is easily computed to give rise to the 
following covariances (QtWt) = g 2 (l — e _7 *)/7 and (Q 2 ) — <? 2 (1 — e~ 2 ~ ft )/2j [30]. As a result, for a constant 
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force the covariance matrix of positions and momenta reads 



(A qt A Pt) = ^-^r^^-^+o^ 

(Ap?) = (38) 

We now confront these terms with the corresponding ones appearing in SVV and SPV. SVV generates the 
following noise terms 

A% = -Q ( h % 

A Ph = e^Qg + Qg (39) 

with covariances given by (Aqf) = g 2 h 2 (l - e^ h ) /2jm 2 = g 2 h 3 /2m 2 + 0(h 4 ), (Aq h Ap h ) = g 2 h(e~ lh/2 - 
e -3 7 h/2)/2 7 m = g 2 h 2 (l - jh)/2m + 0{h 4 ) and (Ap 2 ) = g 2 (l - e~ 2 ~< h )/2~f. Similarly, for SPV we find that 
(Aq 2 ) = g 2 h 3 /4m + 0(h i ), (Aq h Ap h ) = 3 2 /i 2 (l- 7 /i)/2m + 0(/i 4 ) and (Ap 2 ) = g 2 {l- e - 2 ~< h )/2~f. In conclusion, 
for both SVV and SPV positional variance coincides with p8| up to quadratic order, position-momentum 
covariance up to cubic order, and momentum to any order, of the form sought. 

Before proceeding further, we wish to make a couple of comments regarding the derived schemes. The first 
is that in numerical applications, one is usually concerned with the usage of uniform random number generators 
in lieu of the expensive Guassian ones. This is typically possible for the simple Euler scheme [31] in which a 
uniform variate is used to sample a gaussian process up to second moment. In principle, a second order accuracy 
in the propagator would require the usage of a linear combination of uniform variates in order to sample the 
gaussian distribution up to the forth moment. In our treatment, we have used the information up to the second 
moment of the stochastic processes Q t or W t . This circumstance, mirrored by the presence of two random 
variables arise in the Verlet-like class of integrators, allows to employ uniform random number generators in the 
present case. The second comment regards the SVV scheme, whose integration of the deterministic part has 
already appeared in the literature [32], and shown to be equivalent to the one proposed by van Gunsteren and 
Berendsen [23]. In the present approach, however, the noise terms are not imposed to be equal to the matrix 



(38 1, as in previous approaches, but rather emerge spontaneously, ultimately due to the decomposition of the 
underlying vector field. 

The derivation of SVV and SPV is based on the exact integration of eqs. ( 27|28 1. However, one could 



integrate the decomposed dynamics by evaluating the time integrals with expressions analogous to (11 1, but for 
the stochastic case. In particular, by defining the maps 



<7o 

-7/ipo + hF a + W h 



S {-,k Ph+ q0 HF 0+ W h ) ^ 

and the compositions ^ c ,h/2 ° <h,h ° ^c,h/2 or <ki,h/2 ° ^C,h/2 ° ^c,h/2 ° <k(,h/2 another pair of Verlet-like 
algorithms is derived. The updating scheme analogous to SVV, which we name SVVm, then reads 

P* = i + \ h/2 { P ° + \ F{qo) + ^ (1) ( mkT ^ h ) 
h * 

qh = q a -\ — p 



rn 



Ph 



(l-^p^+^Fiq^+^HmkTjh) (41) 



8 



A straightforward calculation shows that the corresponding moments are again accurate to second order with 
the ones of eq.([38|. 

It is easily verified that the deterministic VV and PV propagators are recovered from SVV and SPV (or the 
approximate SVVm) in the limit 7 — ► 0. If morever the Jacobian is phase-space independent, the method is 
called quasi-symplectic (29]. Let us consider, as an example, the stochastic Euler-A method, 



dx h/2 



dx h ( dx h \ / dx \ 1 ( T , h,, n \ ( - G h/2 fe lh l 2 -\ 



and 

'"^-rara .\, + ^)[^-\^f±) C ^) (43, 

so that, for a block diagonal f2, ||J|| = e -' iN i h / 2 _ The same result is derived for the stochastic Euler-B scheme. 
Globally, SVV and SPV have a Jacobian given by||J|| = e~ 3N ' yh . Moreover, this property allows to employ the 
propagators in a Hybrid Monte Carlo scheme, based on an underlying Langevin dynamics. 

The statistical distribution associated to the numerical Langevin flow map is readily derived. Its evolution 
is given by 

f(y, t+h)- f(y, t)= fdx [6(y -x t - Ax) - S(y ~ x t )] f(x t ,t) = f dx { -^S(y - x t )d n ((Ax n )J(x, t)) 

J •* n=l ' 

(44) 

having Taylor expanded the delta function and integrated by parts. Moreover, (•)„ indicates averaging over 
noise. By considering a first order scheme, e.g. the stochastic Euler A, and up to a first order dependence in 
the timestep, we write Ax = Axh + Ax 7 , where Axh is the variation due to Hamiltonian dynamics (at 7 = 0), 
being the vector field arising from the shadow Hamiltonian H of the corresponding Hamiltonian numerical 
scheme. Moreover, Ax 7 arises from the dissipative and random terms. 



By retaining the first two terms of the right hand side of eq. (44 1 we arrive at the following evolution 
equation 

"1 



f(y,t + h) - f(y,t) = -~{Ax H ) v df(y,t) - n { G 



— (Axj) v + mkTd 



f(y,t) (45) 



where we used the fact that, according to standard Langevin analysis, (Ap 2 ) v = mkTh which dominates over 
(Aq 2 ) v ~ h 2 and T is an effective numerical temperature associated to the momentum fluctuations. Moreover, 



in deriving (45 I explicit use of the symplectic character of the Hamiltonian flow has been used. Equation 
([45} represents the numerical Liouville equation coupled to a Fokker-Planck evolution, with stationary solution 
/ oc e-^/ fcf . 

For a second order quasi-symplectic numerical scheme, the information gained from the elementary building 
blocks (Euler A and B) can be used and at equilibrium the Boltzmann distribution is recovered up to order h 2 . 
In this case, the Boltzmann factor contains the shadow Hamiltonian of the underlying scheme at 7 = 0. 

We perform some numerical tests primarily aimed at controlling the quality of the integrators and the 
convergence of the first two moments to the theoretical expectations. An elementary test is provided by the 
one-dimensional stochastic oscillator, with potential U(q) = q 2 /2, integrated in a single timestep scheme with 
the SVV and SPV algorithms. We use h = 0.1, 7 = 0.1 and T = 0.5. The timestep is well below the stability 
limit h = 2, valid for the purely deterministic VV method. The computed momentum and configurational 
distributions, P S i m (p,t) and P S i m (q,t), should converge towards the theoretical forms, P(p, 00) ~ e _/3p I 2 and 
P(q, 00) ~ e~P u ( q \ respectively. The rate of convergence of the normalized histograms is monitored by the 
norms 

E p (T) = a p ma,x\P sm (pi, T) - P(pi, co)\ 

E q (T) = a g max\P sim ( qi ,T)- P(q h oo)\ (46) 

where N is the number of bins and a q and a p are arbitrary factors, as a function of the sampling time window 
T. In Fig. [T]the histograms for the SVV motion are shown, and similar profiles are found for the SPV case. 
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The plot shows that the kinetic and configurational moments are correctly sampled down to the distant tails of 
the distributions. The rate of convergence (Fig. |2j is systematic, and faster for the configurational counterpart. 

We next follow the momentum and position second moments (i.e. the so-called kinetic and configurational 
temperatures) in the high timestep/high friction limit. In this circumstance, it is often reported that Langevin 
integrators produce distinct and systematic departures of kinetic and configurational temperatures from the 
present input, such that equipartion is violated, i.e. (q 2 ) ^ (p 2 ) ^ kT [33 . Here, we employ the information 
from the underlying Hamiltonian propagator to analyze such behavior and showing that, given the quasi- 
symplectic form of the integrator, very accurate results can be obtained as compared to non quasi-symplectic 
ones. The VV update reads 

Qh = (1 - y)<?o + 

Ph = (l-ybo-Ml-^ko (47) 

showing that, to second order, the dynamics arises from a shadow Hamiltonian of the form H — ^p 2 +^(l — \)q 2 . 
In spite of the symmetry of the equations of motion in the q,p variables, the fact that the discretized evolution 
presents a biased, reduced force constant is apparently odd. Moreover, if now the dynamics is equipped with 
the Langevin thermostat, it is understood that equipartition is violated, i.e. positional fluctuations deviate 
quadratically from the kinetic ones [33]. 

However, if we now consider the PV algorithm, constructed via a shift of the updating algorithm by half 
timestep with respect to the VV one, the update reads 

h 2 h 2 
qh = (1 - y)9o + h(l - — )p 

h 2 

Ph = (1 - y )P0 - hq (48) 

showing that, by a simple shift of the "observation" time, the mass of the particle is rescaled by the same factor 
(1 — h 2 /4). Such factor is a manifestation of the symmetry of the original continuous dynamics, although it 
shows up at shifted times as compared to the VV case. If the Langevin thermostat is added, one observes a 
systematic shift of the kinetic temperature from the input one at increasing timesteps. The outcome of the 
present argument is that in order to obtain equilibrated observables, functions of position or velocity separately, 
one should samples these quantities at unequal times, representing optimal sampling points along the trajectory. 

In Fig. [4] we plot the configurational and kinetic temperature of the harmonic oscillator sampled at mid 
((m + h)h, m — 0, 1,2...) and full (mh, m = 0, 1, 2...) time steps. The data show that SVV produces excellent 
averages for both kinetic and configurational temperature at high timesteps, only if these are sampled at unequal 
times. It is important to notice that such averages are extremely robust, in fact they keep close to the theoretical 
value up to timestep h — 2. Raising friction has the effect of shifting configurational temperature at mid-steps 
by about 10%, while kinetic temperature at full-steps weakly deviates from the input value. In conclusion, at 
moderate friction one has a way to sample well equilibrated quantities since the ballistic behavior is well under 
control, while for increasing 7 systematic errors appears in the fluctuations. 

The accuracy of the SVV, SPV and SVVm integrators to sample the dynamical evolution is checked for the 
same stochastic harmonic oscillator for kT — 1 and 7 = 0.1. Given the initial condition q(0) = p(0) = 0, the 
second moments are given by |34] 

(Q 2 (t)} - (q(t)} 2 = 1 + — [-4 + cos%/3i-%/3sinV3t] 

o 

e -t 

(p 2 (t)} - {p{t)) 2 = 1 + — [-4 + cos Vzt + \/3 sin Vst] (49) 



In Fig. [3] we report the numerical relative error of the moments from the theoretical expectation (49 1 for times 
^> h. The convergence of results towards the theoretical values is apparent and follows a h 2 dependence in all 
cases, qualifying the methods as weakly second order accurate. 
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4 Multiple time step dynamics 



In this section, the multiple time step splitting extensions for both the deterministic and stochastic dynamics 
are derived, both with underlying separable Hamiltonian. We consider the decomposition of the Hamiltonian 
into slow and fast components as 

dH(x) = dH s (x) + dH f (x) (50) 

The standard practice in MTS algorithms is to integrate the slow and fast components separately, with timesteps 
h and h/n, where n is a positive integer (n > 1), respectively. A convenient decomposition is given by 



dH s = 



~F s (q) 




dH f = 



-Ff(q) 
p/m 



(51) 



where F s (q) and F s (q) are the slow and fast components of the interatomic forces. The MTS version of the 
VV algorithm is given by 



If f 1 n 

X h = <t>C,h/2° < t>C,h/2n 0( t ) U,h/n°<f>c, h/2n ° <t>C,h/l{ X o) 



(52) 



best known as the RESPA algorithm |8], which is again measure preserving and second order accurate in time. 
An equivalent MTS version of the position Verlet propagator can be constructed along similar lines. 

The stochastic version of the MTS algorithm is obtained by associating the frictional and noise Langevin 
forces to the fast part of the mechanical forces. For example, the MTS version of the propagator generalizing 
the SVV method, reads 



!J C : h/2 ' 



f f 1 71 

^£,h/2n ° <W/n ° W C ,h/2n ° ^C,h/2^o) 



(53) 



We evaluate the performance of the MTS-SVV scheme by considering a one-dimensional model provided 
by a stochastic oscillator with potential energy U(x) = 4.5x 2 + 0.025x 4 . The quartic term of the potential is 
associated to the slow forces. The timesteps are 7r/3 and tt/30 for the slow and fast propagators respectively. 
Friction is set to 1.0 and 0.001 and temperature is set to 0.5. For this choice of the timesteps, the deterministic 
MTS-VV is known to produce a resonance phenomenon |35l IT?] , an effect that limits its applicability in more 
complex situations, such as in presence of intramolecular bond stretching motion. 

Fig. [5] illustrates the convergence towards the expected probabilities of the configurational and momentum 
distributions monitored by the norms (46 1. The results exhibit a systematic convergence for all values of friction 



chosen. Despite its simplicity, the test demonstrates the consistency in the definition of momenta even in a 
MTS type of propagation. 

As a more realistic test, we consider a system made of 256 water molecules and modelled by the flexible 
TIP3P force field (see ref.|9] for details on parameters). Following [36 , the forces are splitted into four different 
levels, although other choices could be made (37]. The first, innermost level handles the fast intramolecular 
bond stretching and angular forces. The second, third and fourth levels integrate the short-range forces arising 
from the Lennard- Jones and direct-space Coulomb interactions (via the Ewald technique [T|) in cut-off regions 
of [0,5], [5,7.5] and [7.5,9.5] A, respectively. Moreover, the fourth level handles the reciprocal space part of 
the Ewald interactions. Friction is set to 7 = 1 ps^ 1 and 7 = 10 3 ps _1 . Temperature is set to T = 300 K. 
The simulations are run with timesteps (h a , nih , n.2h , n^,h ) associated to the four sub-propagators, by setting 
h a = 0.25 fs for the fast propagation and choosing for the integers (m, n 2l 113) the values (1, 1, 1) (STS), (8, 1, 1), 
(8,2,2) and (8,2,4). All simulations are run for 20 ps total time. 

In Fig. [6] the three radial distribution functions for the oxygen-oxygen, oxygen-hydrogen and hydrogen- 
hydrogen atoms are reported for the (1,1,1) and the (8,2,2) runs with 7=1 ps^ 1 . All data converge to the 
correct profiles, i.e. the (1,1,1) run, and equal to a set of independent simulations made by using the purely 
deterministic Nose-Hoover thermostat. However, the (8,2,4) data present distorted profiles for 7=1 ps^ 1 . A 
closer inspection shows that the kinetic temperature of the system is larger than the input value by approximately 
30%, probably due to the large timestep associated to the reciprocal term of electrostatics. However, by simply 
increasing friction to 7 = 10 2 ps^ 1 , the simulation temperature approaches the input value, with relative 
difference being less than 5%, basically removing the distortion in the radial profiles. It should be noticed that 
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with the increased friction we maintain the condition 7/1 < 2 for each level of integration, which ensures that the 
configurational temperature remains close to the kinetic one. In conclusion, some overdamping of the dynamics 
is capable of reducing spurious temperature shifts and configurational bais together with stabilizing the separate 
propagators of the MTS scheme. A separate study should be undertaken to investigate more systematically the 
application of the Langevin MTS approach in the simulation of condensed matter systems. 

5 Conclusions 

The present paper described the integration of Hamiltonian and Langevin dynamical equations guided by a 
symplectic decomposition of the underlying vector field in phase space. In the purely deterministic case, taken 
together with appropriate quadrature formulae, the scheme provided the basis for the Verlet and Leapfrog 
family of numerical propagators, which are second-order accurate and applicable to general deterministic (e.g. 
non-separable Hamiltonian) dynamics. 

In presence of stochastic forces the approach maintains a similar splitting of the underlying Hamiltonian 
dynamics. The vector field route proved convenient for the analytical treatment and we derived numerical 
propagators which are weakly second order accurate. The correctness of the configurational and momentum 
averages and their fluctuations was demonstrated in a series of numerical tests. The deterministic and stochas- 
tic propagators are rapidly extended to multiple time step dynamics, as often employed in the simulation of 
condensed matter systems. In deterministic multiple time step algorithms the presence of resonance phenomena 
limits the stability range of the methodology. By studying a realistic system composed by water, and adopting 
a standard stochastic stabilization with large timesteps, we found that a time-discretization error, manifesting 
itself as a spurious shift in temperature, was reduced by employing a rather large value of the friction coefficient. 

The performances of the proposed schemes to treat more complex stochastic equations, such as dissipative 
particle dynamics [38 , where dissipative and random forces depend in a non trivial way on momenta and 
positions, will be described in a forthcoming paper. 
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Captions 



Figure 1: Histograms P(q,T) and P(p,T) in log-linear scale as a function of z = q 2 and z = p 2 , respectively. 
Results generated with SVV for T — 10 7 h simulation time. The x-axis has been rescaled by an arbitrary factor. 
Squares and circles refer to momentum and position distributions, respectively. 



Figure 2: Norms E p and E q as a function of the sampling time window for the SVV algorithm. The solid line 
refers to E q and the dashed line refers to E p . 



Figure 3: Relative error of second moments from the theoretical expectation (49 1 as a function of the timestep 
sampled for t — 1 and q(0) — p(0) — 0. Sampling time is 10 8 for each timestep chosen. Symbols refer to 
SVV data for positions (open circles) and momenta (filled circles), SPV data for positions (open squares) and 
momenta (filled squares), SVVm data for positions (open diamonds) and momenta (filled diamonds). The dashed 
line corresponds to h 2 and the dot-dashed line to h trends. 



Figure 4: Deviation of configurational and kinetic temperatures from the input temperature for the stochastic 
harmonic oscillator propagated with SVV at kT = 1. Symbols: configurational data sampled at mid-step 
(circles), and full-step (stars), and kinetic data at mid-step (squares) and full-step (triangles). Filled symbols 
refer to 7 = 0.1 and open symbols to 7 = 1.0. 



Figure 5: Convergence towards the expected distribution for the SVV algorithm in the single (thick lines) and 
multiple (thin lines) timestep formulations. Solid and dashed lines refer to configurational and momentum data, 
respectively. 



Figure 6: Radial distribution functions for water for the (1,1,1) (solid curves) and (8,2,2) (dashed curves) 
simulations with the MTS-SVV integrator (see text for details). The curves refer to oxygen-oxygen, oxygen- 
hydrogen and hydrogen- hydrogen profiles. 
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